Density regulation amplifies environmentally induced population fluctuations

Background Density-dependent regulation is ubiquitous in population dynamics, and its potential interaction with environmental stochasticity complicates the characterization of the random component of population dynamics. Yet, this issue has not received attention commensurate with its relevance for descriptive and predictive modeling of population dynamics. Here we use a Bayesian modeling approach to investigate the contribution of density regulation to population variability in stochastic environments. Methods We analytically derive a formula linking the stationary variance of population abundance/density under Gompertz regulation in a stochastic environment with constant variance to the environmental variance and the strength of density feedback, to investigate whether and how density regulation affects the stationary variance. We examine through simulations whether the relationship between stationary variance and density regulation inferred analytically under the Gompertz model carries over to the Ricker model, widely used in population dynamics modeling. Results The analytical decomposition of the stationary variance under stochastic Gompertz dynamics implies higher variability for strongly regulated populations. Simulation results demonstrate that the pattern of increasing population variability with increasing density feedback found under the Gompertz model holds for the Ricker model as well, and is expected to be a general phenomenon with stochastic population models. We also analytically established and empirically validated that the square of the autoregressive parameter of the Gompertz model in AR(1) form represents the proportion of stationary variance due to density dependence. Discussion Our results suggest that neither environmental stochasticity nor density regulation can alone explain the patterns of population variability in stochastic environments, as these two components of temporal variation interact, with a tendency for density regulation to amplify the magnitude of environmentally induced population fluctuations. This finding has far-reaching implications for population viability. It implies that intense intra-specific resource competition increases the risk of environment-driven population collapse at high density, making opportune harvesting a sensible practice for improving the resistance of managed populations such as fish stocks to environmental perturbations. The separation of density-dependent and density-independent processes will help improve population dynamics modeling, while providing a basis for evaluating the relative importance of these two categories of processes that remains a topic of long-standing controversy among ecologists.


INTRODUCTION
Density-dependent regulation is a pervasive feature of population growth processes. It often operates in concert with demographic stochasticity and environmental noise to generate temporal fluctuations in population abundance/density (e.g., May, Hassel & Southwood, 1974;Royama, 1992;Woiwod & Hanski, 1992;Lande, Engen & Saether, 2003;Brook & Bradshaw, 2006). The prevailing form of density regulation is compensation or negative density dependence, which describes a pattern of decreasing population growth rate with increasing population abundance/density, and vice-versa. This phenomenon may result from intraspecific competition for limited resources (e.g., Hansen et al., 1999) or other mechanisms such as predation and diseases that affect net population growth rate in a density-dependent fashion (Hixon, Pacala & Sandin, 2002). Hixon, Pacala & Sandin (2002) highlight the following three salient features of populations undergoing compensatory density regulation. (1) Persistence: the population persists for many generations; (2) boundedness: the population size remains between some positive lower and upper bounds, and (3) return tendency: the population tends to increase below a certain threshold, the equilibrium level, and to decrease when above that threshold. The return tendency is a stabilizing mechanism due to its proclivity to restore populations to their equilibrium levels following a disturbance (Murdoch, 1994;Yodzis, 1995). This feature explains how the harvesting of an abundant population may increase rather than decrease total production in the next generation, and is consequently essential to the concept of sustainable yield in fisheries and wildlife management (Rose et al., 2001;Fowler, 1987).
A relatively less documented form of density dependence is depensation also known as inverse density dependence, positive density dependence or Allee effect (Allee & Bowen, 1932), which refers to a pattern of decreasing population growth with decreasing population density at low densities. This form of density dependence may arise from a variety of mechanisms, including the tendency for predators to kill a fixed number of prey, causing the death rate of the prey population to be higher at low density, and the decrease in birth rate at low population densities due to the difficulty of finding mates. Though rarely detected in practice, the Allee effect is reportedly widespread in nature (Allee & Bowen, 1932;Dennis, 1989;Courchamp, Clutton-Brock & Grenfell, 1999;Kramer et al., 2009). There are at least two reasons for the rare detection of the Allee phenomenon namely, (1) the difficulty of detecting natural populations at low density, and (2) the distortion of statistical analyses by the high variance inherent in small sample sizes (Drake & Kramer, 2011;Kramer et al., 2009). Drake & Kramer (2011 extended the logistic growth model to incorporate the Allee effects. Here we restrict our attention to negative density dependence. There has been since the 1920s a long debate among ecologists on the relative importance of exogenous (environmental) versus endogenous (density-dependent) factors in driving temporal fluctuations in population size/density. This debate reached its pinnacle during 1950s with on the one hand the density-independent and on the other hand, the density-dependent lines of thought led respectively by Andrewartha &Birch (1954) andNicholson (1957). The discovery, by May (1976) that simple discrete-time models of density dependence could generate very complex and potentially chaotic dynamics provided grist to the mills of Nicholson's school, reinforcing the view that observed empirical patterns could be explained without resorting to stochastic factors.
Over recent decades, evaluations of model predictions against observed patterns have resulted in a broad recognition that temporal fluctuations in population abundance/density result from both density-independent and density-dependent factors, potentially interacting in non-trivial ways (Coulson, Rohani & Pascual, 2004). This emerging consensus has shifted the research agenda in population ecology from the simple detection of density-dependence in population time series to the assessment of population dynamical consequences of density regulation. Despite its high relevance for descriptive and predictive population dynamics modeling, the interplay of density-dependent and densityindependent population dynamics processes has not received adequate attention (but see e.g., Fromentin et al., 2001;Ohlberger, Rogers & Stenseth, 2014).
In this study, we analyze, through a combination of analytical derivations and numerical simulations, the interaction between environmental stochasticity and density regulation in driving temporal fluctuations in population abundance. Using the stochastic Gompertz model to describe the population dynamics in a stochastic environment with constant variance, we analytically derive an explicit formula linking the stationary variance of population abundance/density to the environmental variance and the strength of density feedback. We derive a formula quantifying the contribution of density regulation to population variance in stationary phase. We conduct a simulation study to corroborate empirically the analytically established relationship between density feedback and stationary variance and to check whether the same pattern of association holds for the Ricker model (Ricker, 1954), which is also widely used for population dynamics modeling, particularly in fisheries.

Model specification and analytical derivations
Let Y 0 and Y t denote respectively the initial population size (population size at time 0) and the population size at time t (t ≥ 1). We assume a discrete-time stochastic Gompertz model for the population dynamics (Reddingius, 1971;Royama, 1981;Sibly et al., 2005;Dennis et al., 2006;Mutshinda & O'Hara, 2011;Mutshinda, O'Hara & Woiwod, 2009;Mutshinda, O'Hara & Woiwod, 2011;Mutshinda et al., 2019) so that the population size at time t ≥ 1 is given by where r > 0 is the intrinsic growth rate (i.e., e r is the multiplicative population growth rate in the absence of density dependence), whereas the realized growth rate Y t /Y t −1 decreases with increasing population size Y t −1 , and k > 0 is the natural logarithm of the carrying capacity or equilibrium population size/density denoted by K . The error term ε t , assumed to be Gaussian with mean zero and variance σ 2 t , typically integrates process stochasticity (demographic and environmental stochasticity) and observation error. The Gompertz model as presented here is phenomenological and only describes population level changes without explicitly accounting for individual level processes or age structure, in contrast with mechanistic models that translate individual parameters into population level consequences.
For the purpose of the present study, we rely on simulated data and so, we assume that the population size is observed without error. In addition, we know the data-generating model when in practice the error term typically involves the effect of model misspecification as well. For small populations, demographic stochasticity may be important, and its inverse scaling with the population size introduces time-dependence in the random error term σ 2 t (e.g., Saether et al., 2000). Since our focus is on the population dynamics in stationary phase, we assume that the population size is large enough for demographic stochasticity to be unimportant. Therefore, we consider the process error to be entirely due to environmental fluctuations in a stochastic environment with constant variance, so that σ 2 t = σ 2 with the error terms ε t , t ≥ 1 being serially independent. On the natural logarithmic scale with y t = ln(Y t ), model (1) becomes Letting β = 1 − rk −1 , we can re-write model (2) as which is a first-order autoregressive [AR(1)] model with autoregressive coefficient β. From the expression of β in terms of r and k, it follows that k = r(1 − β) −1 , which establishes a one-to-one correspondence between the parameters of the Gompertz model in AR(1) form (Eq. 3) and those of the standard Gompertz model (Eq. 1). However, the AR(1) form is more convenient for analyzing the model's dynamic behavior as we do next.
We start by analyzing the dynamic behavior of the deterministic Gompertz model in AR(1) form (i.e., Eq. 3 without the error term ε t ). We discuss the existence of an equilibrium point and the pace at which the population size/density approaches the equilibrium. We then consider the stochastic version of the model (Eq. 3) for which the equilibrium is a non-degenerate distribution (the stationary distribution) rather than being a single point as in the deterministic model. We analyze the temporal fluctuations of the population size around the mean value of its stationary distribution.

Dynamic behavior of the deterministic Gompertz model
The deterministic part of the Gompertz model in AR(1) form is given by the first-order difference equation If |β| < 1, the system will eventually converge to an equilibrium state y ∞ at which the population size remains constant, meaning that y ∞ = r + βy ∞ . Solving this equation for y ∞ yields the following expression for the non-trivial equilibrium We can use the relation r = y ∞ (1 − β) resulting from (5) to describe the population dynamics in terms of y ∞ as y 1 = y ∞ +β y 0 − y ∞ , y 2 = y ∞ +β 2 (y 0 −y ∞ ), and by induction we obtain the following closed-form solution to the first-order difference Eq. (4).
Consequently, y t will converge to y ∞ provided that |β| < 1, with faster convergence when |β| is closer to zero, with y t approaching the equilibrium state y ∞ either monotonically from above or from below when 0 < β < 1 or through damped oscillations when −1 < β < 0. Following from Eq. (6), y t relates to y t −1 as y t = y ∞ + β y t −1 − y ∞ , which can be re-arranged as y t = (1 − β)y ∞ + βy t −1 . The relationship y t = (1 − β)y ∞ + βy t −1 indicates that y t is a linear combination of y ∞ and y t −1 with the autoregressive coefficient β providing a direct measure of the dependence of y t on y t −1 . As a result, the autoregressive parameter of the Gompertz model in AR(1) form is often referred to as the strength of density dependence (e.g., Hampton et al., 2013;Ponciano, Taper & Dennis, 2018;Messmer, 2019;Peeters et al., 2022). We also adopt this terminology.
Deterministic models of population dynamics draw on the assumption that vital rates are constant over time, so that a single set of input values uniquely determines the output value. However, an obvious feature of the real world is that the environment varies continually. Environmental stochasticity refers to the variability in demographic rates caused by random variations in environmental conditions, whereas sampling variations in independent outcomes of demographic events (births, death and dispersal) among individuals in a finite population produces a different kind of fluctuations in population dynamics known as demographic stochasticity. While environmental stochasticity affects all populations irrespective of size, demographic stochasticity is generally only relevant in small populations due to its inverse scaling with the population size. We refer to Lande, Engen & Saether (2003) for details on statistical methods for estimating demographic and environmental stochasticity from empirical data. We next consider the dynamic behavior of the stochastic Gompertz model and investigate the interaction between density-dependent effects and environmental noise in driving temporal fluctuations in population abundance/density.

Dynamic behavior of the stochastic Gompertz model
The stochastic Gompertz model in AR(1) form is given by Eq. (3). If |β| < 1 the population size will eventually reach a stationary distribution (or equilibrium distribution), which is the stochastic version of an equilibrium in the deterministic model. In stationary phase, the population growth rate is zero on average. Since the additive error ε t are assumed to be normally distributed, the stationary distribution is also normal with long-run mean value µ ∞ and variance v ∞ determined by their time-invariance property. Taking expectations of both sides of Eq. (3) and applying the time-invariance property of the long-run mean yields the equation µ ∞ = r + βµ ∞ whose solution for µ ∞ is Alternatively, taking variances of both sides of (3) and applying the time-invariance property of the long-term variance yields the equation Equations (7) and (8) indicate that both the mean of the stationary distribution, which is identical to the deterministic stable equilibrium value y ∞ = k (the log-carrying capacity), and the stationary variance (i.e., the variance of population time series in stationary phase) depend on the strength, β, of density regulation. In addition, the stationary variance v ∞ is proportional to the environmental variance, σ 2 , and relates to the density feedback in such a way that stronger density regulation (i.e., |β| values close to 1) induces higher variability.
It is worth emphasizing that stationarity is a key requirement for analyzing density dependence in population time series, and Eqs. (7) and (8) only make sense for stationary time series. Therefore, the sensible range for the autoregressive parameter β is the interval (−1,1). However, since we are interested in the stationary variance which depends on β only through β 2 , we restrict attention on β values in the unit interval 0 ≤ β < 1. From our model assumptions that r and k are both positive and the expression β = 1 − rk −1 connecting the autoregressive coefficient of the Gompertz model in AR(1) form to the parameters r and k of the standard Gompertz model, it follows that 0 ≤ β < 1 corresponds to r ≤ k.
When β = 0, which corresponds to the case r = k in the standard Gompertz model formulation, the population trajectory follows a Gaussian white noise process shifted at r. This process is stationary with constant variance σ 2 . When β = 1, which arises as a limiting case when k tends to infinity in the standard Gompertz model formulation, the dynamics are density independent, and the population trajectory is a random walk process with drift r. This process is not stationary since Var y t = t σ 2 is unbounded and E y t = y 0 + rt is only constant when r = 0 (i.e., for the random walk without drift), while a key condition for (weak) stationarity is that the mean and variance of the time series are constant, the other requirement being that for any integer h, Cov(y t ,y t +h ) depends only on h and not on t .
For density-dependent dynamics where β is statistically different from zero and |β| < 1, the population fluctuates around the stationary mean µ ∞ with variance σ 2 /(1−β 2 ), where σ 2 represents the environmental variance. Therefore, the portion of the stationary variance due to density regulation, herein denoted by σ 2 dd , is given by the difference between the stationary variance σ 2 /(1 −β 2 ) under density-regulation and its counterpart σ 2 in the absence of density feedback on population growth rate. That is, Consequently, the proportion of the stationary variance v ∞ due to density regulation is Interestingly, the proportion of the stationary variance v ∞ due to density regulation is simply the square of the density-dependence effect β. This provides a biological interpretation of the autoregressive coefficient of the AR(1) model when applied to population dynamics, besides being a measure of the strength of density dependence. Equation (9) provides a tool for separating the stationary variance of population time series into a density-independent component and a density dependent counterpart.
A key motivation for using the Gompertz model is that it can be written as an AR(1) model, which is simple with known statistical properties and established inferential procedures (Mutshinda & O'Hara, 2010). Because the AR(1) model is a linear model with Gaussian errors, maximum likelihood estimates of r, β , and σ 2 , which are identical to least square estimates, can be obtained by performing a linear regression of y t on y t −1 (t = 1,2,3,...) using standard statistical packages. However, confidence intervals from standard statistical packages are not valid because the y t s are not independent due to the autoregressive structure in the model. Dennis & Taper (1994) discuss bootstrap and jackknife procedures for constructing confidence intervals of the AR parameters. In a Bayesian framework, priors distributions can be defined to appropriately constrain the model parameters to more sensible range of values, thereby reducing posterior uncertainty while guaranteeing parameter identifiability (Mutshinda, 2009a;Mwanza, 2010).
A common objection to the Gompertz model is that the growth rate depends, if applicable, only logarithmically on the population density allegedly inducing weaker density feedback than the closely related Ricker model (Ricker, 1954) where the growth rate at time t depends on Y t −1 , the actual population size (Dennis & Taper, 1994). The Stochastic Ricker model is given by where r is the intrinsic growth rate and K is the carrying capacity or equilibrium population density on the scale of untransformed population size, and ε t are zero-mean random shocks to the population growth rate assumed to be normally distributed and serially independent. The Ricker model is widely used for population dynamics modeling, particularly in fisheries.
In the next section, we describe a simulation study designed to corroborate empirically the pattern of increasing stationary variance with increasing strength of density regulation under the Gompertz model and to investigate whether this pattern carries over to the Ricker model.

Data simulation
Given numerical values of the parameters r,k and σ , an initial log population size y 0 , and a routine for generating standard normal random variables, one can recursively generate trajectories y 0 , y 1 , y 2 , . . . y n from the Gompertz model with initial population density set to the log-carrying capacity. We simulated replicated population time series over 300 time steps and discarded the first 200 observations to ensure that the last 100 data points come from the stationary distribution. We tuned the parameters of the data-simulating model to mimic different levels of density regulation and different levels of environmental noise. Since the parameter β only depends on the positive parameters r and k through the ratio r/k , one can generate data with different levels of density-regulation (different β values) by fixing one of the two parameters, typically the carrying capacity, and varying the other (e.g., Greenwell & Ng, 1984). Fixing the carrying capacity to 1 amounts to expressing population density in units of the carrying capacity. In our simulations, we fix the log-carrying capacity to 1, so that β = (1 − r) for 0 < r < 1. We simulated data under four different levels of environmental noise, with environmental variance set to 0.10, 0.15, 0.20 and 0.25, and three different values of r namely, 0.8, 0.6, and 0.4, corresponding to β values 0.2, 0.4, and 0.6, respectively.
We graphically checked whether the empirical stationary variance v ∞ increased with increasing the strength of density feedback as implied by the analytical result ( (8)). We computed the contribution σ 2 dd of density regulation to the empirical stationary variance v ∞ at different levels of density regulation as the difference between the empirical stationary variance of simulated population trajectories and the assumed environmental variance.
In practice, the model parameters are unknown, and one has to rely on estimates obtained from the model fitting to data. We examined the relationship between the strength of density regulation and the proportion of stationary variance due to density regulation using environmental variance estimates over 100 replicated population trajectories at each of the three levels of density regulation under consideration.
A question that springs to mind is whether the pattern of increasing stationary variance with increasing strength of density feedback inferred under the Gompertz model holds for the Ricker model as well. Since under the Ricker model, unlike the Gompertz model, we do not have an expression separating the stationary variance into contributions from environmental stochasticity and density regulation, we rely on simulations to tackle this question. We fit, with a Bayesian approach (Gelman et al., 2013;Mutshinda et al., 2022), the stochastic Gompertz and stochastic Ricker models to data replicates simulated from the stochastic Gompertz model with different levels of density regulation. For the Gompertz model, we independently assigned a Gamma(1,1) prior on the intrinsic growth rate r, a standard normal prior independently on the autoregressive coefficient β and an InvGamma(0.1,0.1) prior on the environmental variance σ 2 . For the Ricker model, we independently assigned a Gamma(1,1) prior on the on the intrinsic growth rate r and Gamma(0.1,0.1) priors on the carrying capacity K and the environmental variance σ 2 . We used Markov chain Monte Carlo methods (Gilks, Richardson & Spiegelhalter, 1996) to simulate, via the Bayesian freeware OpenBUGS (Thomas et al., 2006;Mutshinda, 2009b), from the joint posterior distributions. We primarily ran the models in OpenBUGS to assess their convergence both informally through visual inspection of traceplots and autocorrelation plots in OpenBUGS, and formally by looking at the Gelman, Rubin statistics via the R package CODA (Plummer et al., 2006), and found that at least 2,000 iterations were required for convergence. Therefore, we used a burn-in period of 4,000 iterations for both models. We estimated the contribution σ 2 dd of density dependence to the stationary variance v ∞ by the difference between the empirical stationary variance of simulated population time series and the posterior estimate of the environmental variance, and subsequently derived the proportion ϕ dd of stationary variance due to density regulation as ϕ dd = σ 2 dd /v ∞ .

RESULTS
The stationary variance of simulated population trajectories increased monotonically with the strength of density regulation (Fig. 1). This result indicates that in a given environment, strongly regulated populations will exhibit larger temporal fluctuations than weakly regulated ones. In addition, the proportion 1 − σ 2 /v ∞ of stationary variance due to density regulation increased monotonically with the strength of density dependence in simulated population time series (Fig. S1 in Online Supplemental Material). As noted earlier, the environmental variance σ 2 is unknown in practice and needs to be estimated from data. Our Bayesian model fitting was effective at retrieving the environmental variance in population time series under both the stochastic Gompertz and stochastic Ricker models (Fig. S2). This allowed us to evaluate the proportion ϕ dd = σ 2 dd /v ∞ of stationary variance due to density regulation through Eq. (9). Under either model, the proportion of stationary variance due to density regulation increased monotonically with the strength of density feedback in the data (Fig. 2).
In addition, the inferred proportion of stationary variance due to density dependence was approximately equal to the square of the autoregressive parameter β over simulated data replicates for the stochastic Gompertz model (Fig. 3), consistent with our analytical derivation in Eq. (10). This result provides another biological meaning to the autoregressive coefficient of the AR(1) model when applied to population dynamics.

DISCUSSION
In this study, we combined analytical derivations and numerical simulations to analyze the interplay between environmental noise and density regulation in driving temporal fluctuations in population abundance/density. We derived a formula (Eq. 8) relating the stationary variance of the population abundance/density under Gompertz-type density regulation in a stochastic environment with constant variance to the environmental variance and the strength of density dependence, implying that density regulation amplifies the magnitude of environmentally induced population fluctuations. We worked out a formula separating the stationary variance in population abundance/density into its density-independent and density dependent components. An important result emerging from this variance decomposition is that the square of the autoregressive coefficient of the Gompertz model in AR(1) form represents the proportion of stationary variance due to density regulation (Eq. 10).
Simulation results substantiated empirically the analytically established pattern of increasing stationary variance with increasing strength of density regulation under the Gompertz model (Fig. 1). The Bayesian model fitting was effective at retrieving the environmental component of the stationary variance under both the stochastic Gompertz and stochastic Ricker models (Fig. S1). This allowed us to estimate the portion of stationary variance due to density-regulation by the difference between the stationary variance of simulated population trajectories and the posterior estimate of the environmental variance under either model.
We expect the pattern of increasing stationary variance with the strength of density feedback established analytically and/or empirically under the stochastic Gompertz and Ricker models (Fig. 2) implying higher temporal variability for strongly regulated populations to be a general phenomenon with population dynamical models. This is because the linear relationships in the Gompertz and Ricker models can be considered as Taylor approximations near equilibrium of more complex density-dependent growth functions (Dennis & Constantino, 1988).
Simulation results also substantiated empirically the analytically established result in Eq. (10) that the square of the autoregressive coefficient of the Gompertz model in AR(1) represents the proportion of stationary variance due to density regulation (Fig. 3). This finding provides another biological meaning to the autoregressive coefficient β of the Gompertz model in AR(1) form: Besides being a measure of the strength of density dependence, β represents the proportion of stationary variance due to density dependence.
Our main finding that density regulation amplifies environmentally induced population fluctuations has important implications for population viability. It suggests that intense intra-specific resource competition increases the risk of environment-driven population collapse at high density, lending support to opportune harvesting as a means to improve the resistance of managed populations such as fish stocks to environmental perturbations.
Overall, our analytical and empirical analyses demonstrate that the impact of density regulation on population dynamics involves a deterministic component affecting the mean population abundance/density and a stochastic component affecting the process variance. Therefore, the de facto decomposition of process variance into environmental stochasticity and demographic stochasticity disregards the contribution of density dependence to population variance in stochastic environments, which may be substantial.

CONCLUSIONS
We analytically established and empirically verified that environmental noise interacts with density feedback in convoluted ways, causing density-regulated populations to undergo stronger fluctuations than expected under the sole influence of environmental stochasticity. The separation of exogenous (environmental) and endogenous (density-dependent) components of population variability is essential to evaluating their relative importance, which remains a topic of continuous debate among ecologists (Coulson, Rohani & Pascual, 2004).
Our analyses demonstrate that one can effectively extract the environmental component of temporal variability in population dynamics variability by fitting ecological models to population time series, and the Bayesian approach adopted here has proven fruitful to this end.
In order to decompose the stationary variance of population size in its densityindependent component and density-dependent components, we made the following two simplifying assumptions: (1) Demographic stochasticity is unimportant, and (2) population sizes are recorded without error. While the irrelevance of demographic stochasticity in large populations is a sensible assumption due to the inverse scaling of demographic stochasticity with the population size (Lande, Engen & Saether, 2003), real-world observational population time series are typically fraught with observation or sampling errors. Inadequate handling of sampling error can induce wrong inferences about population processes, including spurious density dependence detection (Dennis & Taper, 1994;Freckleton et al., 2006). Density-dependent population dynamics models that integrate process variation and observation error are required to analyze the interaction of environmental stochasticity and density regulation in real-world populations. These kinds of models can be conveniently developed and fitted in the Bayesian state-space modelling framework (e.g., De Valpine & Hastings, 2002;Buckland et al., 2004;Dennis et al., 2006).

ADDITIONAL INFORMATION AND DECLARATIONS Funding
This work was supported by the Simons Collaboration on Computational Biogeochemical Modeling of Marine Ecosystems (CBIOMES) (Grant ID: 549935). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.